Closes #122, #123

This PR implements a variant of the approach hashed out in the above issues (with some significant variations), designed to

• Gracefully deal with degeneracies, including topologically non-trivial ones such as Dirac points
• Allow O(1) access of simplices in a bandstructure

I'm rather happy with the design including here. It is both cleaner, more correct and (sometimes) faster than what we have in master regarding degeneracy resolution in bandstructure.

Degeneracies

For context, master deals with d-order degeneracies by keeping d copies of the degenerate eigenvalues, and co-diagonalizing their eigenvectors with some kind of a velocity operator (be it using DualNumbers or finite differences). The codiagonalizer thing was tricky and rather inelegant (somewhat "heuristic" if you're being kind, an ugly hack otherwise). The approach completely failed to give properly connected bands in the presence of Dirac points and/or full band degeneracies (with or without a non-trivial non-Abelian connection). The essential shortcoming is that the codiagonalization approach was an attempt to generalize the 1D UntangledBands algorithm in MathQ to a higher-dimensional setting. Which, it turns out, is the wrong way to go about this.

The design in this PR takes a different approach that is free of any sort of codiagonalization step. In essence, it keeps degeneracies as single vertices in the bandstructure. Degenerate vertices are associated to a single energy but to a set of eigenstates that form the basis of the degenerate subspace. No attempt is made to choose any specific basis within degenerate subspaces. We stick to what the diagonalizer library produces. Actually, all these subspace basis states are in fact views into the original arrays produced by the diagonalization. Since v1.5 this is allocation-free and very efficient.

Extracting the connectivity of a network of eigenpairs that can exhibit degeneracies is a very interesting exercise. It boils down to just two rules:

• Connect adjacent vertices with the highest projection, where projection is understood as full projection between potentially degenerate subspaces
• The number of connections spawned by a subspace at base momentum k onto the different target subspaces at another base momentum k´ is not simply one. It is the difference in degeneracies. In other words, the degeneracies of all target subspaces should add up to the minimum number equal or greater to the degeneracy of the source.

These two rules produce correct connectivities with any degeneracy configuration that I have tested, including Dirac points, band crossings, Riemann folds and degenerate subbands, as long as the mesh is sufficiently fine to resolve connections through adjacent subspace overlaps.

A consequence of these rules is that any two subbands sharing a degenerate vertex will actually be the same subband. As a result, band folding does no longer split bands into different subbands. Nether do Dirac points, i.e. pi-bands in graphene form a single band that is double-valued on the Brillouin zone.

Preparing for optimized O(1) bandstructure access

An interdependent change is included in this PR. The mesh command is renamed once more (it produced conflicts with Makie, in particular). It is now called cuboid, which is a special kind of (lazy) CuboidMesh, different from BandMesh, which is more similar to the old Mesh. Cuboid meshes are marching-tetrahedra meshes over hyper-rectangles (or cuboids) in D-dimensions. All k-point or parameter sampling domains for bandstructures are of this specific kind now. We can still use all the mapping machinery to convert the cuboid to an arbitrary mapped domain, but restricting the base mesh (as it is now called) to a CuboidMesh opens the door to two future improvements:

• Fast simplex access
• Iterative mesh refinement

More on that when the time comes. But for now, the user has three things to learn about this change

• Write cuboid instead of mesh
• Write subticks instead of points for the cuboid kwarg (points could be mistaken for the total number of vertices in the mesh, hence the change)
• You can now use an extended syntax cuboid((x1, x2, ... xN), (y1, y2, ... yM)....; kw...) where xi, yi etc are ticks on the x, y axes of the cuboid, such that a vertex will be placed at (xi, yj, zk,...) for any i, j, k combination. The subticks are additional equidistant divisions within each pair of adjacent ticks. In other words, we can still do

mesh = cuboid((-pi,pi), subticks = 13)

but now we can also do

mesh = cuboid((-pi, -pi/3, pi/3, pi), subticks = 13)
mesh = cuboid((-pi, -pi/3, pi/3, pi), subticks = (5, 15, 5)

Performance and correctness

The performance of bandstructure relative to master is in general rather similar. This PR actually does quite a bit more work than master, although it does so quite a bit more efficiently. The number of allocations, in particular, is drastically reduced. For comparison we have, with h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) |> unitcell(2)

master:

julia> @btime bs = bandstructure(h, Quantica.mesh((-pi, pi), (-pi, pi), points = 7))
  1.687 ms (6457 allocations: 1.23 MiB)
Bandstructure{3}: collection of 2D bands in a 3D space
  Bands        : 5
  Element type : scalar (Complex{Float64})

(Note the missing simplices around Dirac points)

PR :

julia> @btime bandstructure(h, cuboid((-pi, pi), (-pi, pi), subticks = 7))
  1.337 ms (912 allocations: 1.02 MiB)
Bandstructure{2}: bands of a 2D Hamiltonian
  Bands         : 1
  Element type  : scalar (Float64)
  Vertices      : 310
  Edges         : 924
  Simplices     : 588

(This PR correctly captures the Dirac points)

Other larger systems can actually show a significant performance increase, like this system with subband degeneracies from a project we're working on.

master (timing second run)

julia> @time bs = bandstructure(hN(), Quantica.mesh((-pi, pi), (-pi, pi), points = 80))
  3.201837 seconds (3.54 M allocations: 1.553 GiB, 16.62% gc time)
Bandstructure{3}: collection of 2D bands in a 3D space
  Bands        : 1
  Element type : scalar (Complex{Float64})

(Band degeneracies makes master choke connecting vertices here)

PR (timing second run)

julia> @time bs = bandstructure(hN(), cuboid((-pi, pi), (-pi, pi), subticks = 80))
  0.727641 seconds (104.77 k allocations: 178.502 MiB, 8.72% gc time)
Bandstructure{2}: bands of a 2D Hamiltonian
  Bands         : 1
  Element type  : scalar (Float64)
  Vertices      : 25592
  Edges         : 75540
  Simplices     : 49732

Band degeneracies are handled correctly here. The remaining holes are either inevitable (due to vertices of Riemann folds or Dirac cones that do not fall on a base mesh vertex) or the result of insufficient base mesh resolution at weakly avoided band crossings. The former are essential connectivity frustrations that (I think) can only be resolved by refining the mesh to get the dislocations to fall on a base mesh vertex.

Other details

• Both the Makie and VegaLite interfaces have been updated to plot the new BandMeshes
• All simplex vertices and their corresponding subspace bases are collected into vectors (the new Simplices type) that can be used to efficiently compute density of states and Green's functions in the future.